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Abstract. We report an extension of the smoothed profile method (SPM) [Y. Nakayama, K. Kim, and 
R. Yamamoto, Eur. Phys. J. E 26, 361 (2008)], a direct numerical simulation method for calculating the 
complex modulus of the dispersion of particles, in which we introduce a temporally oscillatory external 
force into the system. The validity of the method was examined by evaluating the storage G'{uj) and 
loss G"(oj) moduli of a system composed of identical spherical particles dispersed in an incompressible 
Newtonian host fluid at volume fractions of # = 0, 0.41, 0.46, and 0.51. The moduli were evaluated at 
several frequencies of shear flow; the shear flow used here has a zigzag profile, as is consistent with the 
usual periodic boundary conditions. The simulation results were compared with several experiments for 
colloidal dispersions of spherical particles. 



PACS. 82.70.-y Disperse systems; complex fluids 
83. 60. Be Linear viscoelasticity 

1 Introduction 

The viscoelastic properties of the dispersion of solid par- 
ticles arc of particular importance in several scientific, 
engineering, and industrial fields. These properties are 

Send offprint requests to: 
a tiwashitQutk. edu 
b ryoichi@cheme.kyoto-u.ac.jp 
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strongly affected not only by direct interactions between 
particles, but also by thermal fluctuations of the system 
and the hydrodynamic interactions (HI) acting on dis- 
persed particles, mediated by the surrounding fluid [112] , 
The viscoelasticity of materials are most commonly char- 
acterized by the complex modulus, which represents the 
mechanical response of materials to temporally oscillat- 
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ing small shear deformations of frequency w. The complex 
modulus consists of elastic and viscous components (i.e., 
the storage modulus G'(uj) and the loss modulus G"(oj)). 
In general, materials tend to lose elastic energy by dissi- 
pation when the frequency of external deformation is low, 
but can store elastic energy when the frequency is high. 

Understanding viscoelastic behavior has been the sub- 
ject of both fundamental and technological interest. Es- 
pecially, the linear and non-linear responses of colloidal 
systems, including colloidal glasses and aggregated col- 
loidal suspensions and gels, have attracted much atten- 
tion in recent year |3I4I28| . While extensive experimental 
studies have been carried out to determine the viscoelas- 
tic properties of particle dispersions over a wide range 
of volume fractions, from dilute ~ 0) to very dense 

> 0.5), theoretical studies encounter a fundamental 
problem when the density of dispersion is high, because 
of many-body effects and the long-range nature of the HI 
among dispersed particles. Developing analytical theories 
becomes even more difficult when the host fluids are com- 
plex, such as electrolytes or viscoelastic media. Computer 
simulations can thus be very powerful tools in theoretical 
investigations of dense dispersions in general. 

The widely used simulation technique is Stokesian dy- 
namics (SD) method [5], which is based on the Stokes 
approximation (Re — s> 0), involves near- field lubrication 
forces and far-field many-body His. The SD method has 
succeeded in simulating the motions of particle disper- 
sions at steady shear. However the viscoelasticity of parti- 
cle dispersions at oscillatory shear has not been examined. 
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Furthermore, the SD method is valid only for simulating 
the particles in a Newtonian fluid and is not applicable to 
the motions of particles in complex fluids, such as charged 
particles in electrolytes or viscoelastic media. 

Various numerical methods have recently been pro- 
posed for simulating dense dispersions, including some 
cases where the host fluids are complex. In order to accu- 
rately track the motions of host fluids as well as the mo- 
tions of dispersed particles, a variety of numerical methods 
have been developed. One of these methods, which we call 
the direct numerical simulation (DNS) method for particle 
dispersions, involves solving the Navier- Stokes (NS) equa- 
tion for host fluids in a manner consistent with boundary 
conditions defined according to the particle positions and 
motions. Aside from several successful implementations of 
DNS methods for particle dispersions 6 7 8 9 1101111121131b 
other methods have been developed that do not rely on 
the Navier-Stokes equation to resolve fluid motions. The 
coupling methods of particle dynamics with the lattice- 
Boltzmann (LB) method [15] or stochastic rotational dy- 
namics (SRD) |16ll7j are the most popular alternatives 
to DNS, because LB and SRD are believed to be more 
computationally efficient. 

Problems studied here have very small Reynold num- 
ber and the full NS equations for the host fluid are not 
necessarily solved to simulate the motions of particles in a 
Newtonian fluid. It is sufficient to use the methods based 
on the Stokes approximation. However there are several 
advantages for solving the NS equations for the host fluid. 
Those techniques mentioned above enable us to simulate 
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the short-time motions of particle dispersions, where the 2 Simulation method 

coupling between the particle motion and the fluid motion 

Let us consider a dispersion composed of identical spher- 

rcmains strongly, even if Reynolds number is very small. 

ical particles of radius a in a Newtonian host fluid sub- 
One of the typical behavior is the power-law decay in the 

jected to oscillatory shear. The host fluid is described by 

velocity correlations of particles, which is known as "long- 

the velocity field v(r,t) and the pressure field p(r, t). The 

time tail" . In addition, Those methods have succeeded in 

ith dispersed particle is described by {Ri, Vi, fli), where 

simulating several complex phenomena, such as the struc- 

Ri is the position of the particle, Vi is the translational 

ture formation in colloidal gels where the pressure filed of 

velocity, and fli is the rotational velocity. The coupling 

the fluid plays an important role |18j and electrophoresis 

scheme between the fluid motions and the particle mo- 
ot charged colloidal dispersions |19l20j . 

tions is based on the smoothed profile method (SPM), 
Although extensions of those methods have already which introduccs a P articlc dcnsit y ficld ^,t) ^ [0 : 1] 
been developed for steady shear flow with DNS ESj, LB on thc cntirc ficld ' ^ = 1 for thc P articlc d °mains and 
j2l l 22 l 23] , and SRD [M] , there exists no successful attempt ^ = for the fluid d °mains. These domains are separated 
for calculating the complex moduli of thc dispersions with h ? a thin intcrfa cial domain of thickness £. Thc SPM is an 
an imposed oscillatory shear flow. The purpose of thc efficient mcthod to resolve the HIs between thc fluid and 
present paper is to report our successful reproductions P article motions ; ^ details are given in [9J10112j. 
of experimentally observed behaviors of the storage G' {u>) The time evolution of the host fluid is governed by the 
and loss G"(u) moduli in temporally oscillating shear flow. Navi cr-Stokcs equation 

Pf (d t v + vVv)=V-<t + Pf 4>f p + Pf f shear (1) 

To check thc validity of the method, we first apply the 
mcthod to a Newtonian host fluid of viscosity 77, which with the incomprcssibility condition V • v = 0, where pf 
should exhibit a purely viscous response, G'{tu) = and is the density of the fluid, rj is the shear viscosity, the 
G"(uj) = r]U). We then apply the method to dense disper- stress tensor a = —pi + r]{Vv + (Vt>) T }, and f shear (r^ t) 
sions composed of identical spherical particles in a Ncwto- is an external force field that is introduced to enforce an 
nian host fluid. Since several experimental measurements oscillatory shear flow on the entire system. The body force 
of the frequency dependence of G'{uS) and G"(uj) for dense (f>f p is introducd to ensures the rigidity of particles and 
colloidal dispersions have already been reported [58 29 30 31i,thc appropriate boundary condition at the fluid/particle 
we compare our numerical results with those experiments, interface |9HOI12j . 

to examine the validity of our method for dense particle The time evolution of the ?th particle with mass Mj 
dispersions. and moment of inertia Ii is governed by Newton's equa- 
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tions of motion: imental time scales. On the other hand, in the simula- 

„ r , tions based on the DNS approach, the acceleration term 

MiVi = ff + ff + G\ , Ri = Vi, (2) 

strongly affects the apparent stress of the dispersion. 

Ii ■n i ^N l H + Gf 2 , (3) 

There are two key points for calculating the apparent 

where F, H and N, H are the hydrodynamic forces and torques, , aOT) r-,-, j- • -\ i , i i , ,i i 
1 1 J J H 'stress a pp of the dispersion: i) how to calculate the accel- 

respectively, exerted by the host fluid on the particle. Ff eration term ^ ^ and ^ how tQ introduce the external 
is the repulsive force that prevents the particles from over- forcc /sW . The firgt pQint ig straightforward . For a sim _ 

lapping, and a truncated Lennard- Jones potential, V(ra) — i , i ,i i ,• , u • , 

I ^ °' 1 uiation time step h. tne acceleration term can be simply 

4[( ( rM J ) 36 -W^) 18 +l/4]forr iJ < 2 1 / ls aovV(r ij ) - 0, calculated as 1/ V J dxx(p t v(x,t + h) - p t v(x,t))/h. Next 
is adopted in this work. Here, a = 2a and r« = |^ - H, |. thc external forcc /8 / l ear is introduccd as a body forcC5 

G/ and Gy arc the random forcc and torque, respectively, , c ,i c n -n i. i ^ c u j-i 

II H ' ^ J ' to enforce the following oscillatory velocity held over the 



due to thermal fluctuations. These random fluctuations 
are assumed to be (Gf ) = (Gf } = 0, (Gf (t) ■ Gf (0)) = 
3a v 5{t) and (Gf (t) • Gf (0)) = 3a n S(t), where (} denotes 
time averaging, and ay and a a are numerical parameters 
to control the particle temperature T. The procedure for 
determining the temperature is described in [25] . 

The apparent stress a app of the dispersion is written 



entire system, 

7(t)(-!/ - L y /2), (-L y /2 < y < -Ly/A) 

<{V) = I 7 (i)y, {-Ly/A <y< Ly/A) (5) 

j(t)(-y + L y /2) (L y /A<y<L y /2) 
j(t) = 70 cos(wt) (6) 



as 



where y denotes distance in thc velocity-gradient direc- 
tion, and 7(i) is the oscillatory shear rate, with amplitude 
* J V J dt ^y and frequency uj. Here, the flow is imposed in the x 

where p t = (1 - <fi)pf + 4>p p , p p is the density of the par- direction and L y is the length of the system in the y di- 



r a PP = 1 J dxxptf sHear _ W dxx ^. {ptv) (4) 



tides, and V is the volume of the system. The derivation rection. This zigzag velocity profile was first used to simu- 

of Eq. g]) was reported in late dispersions in a steady shear |25l26j . Note that f shear 

Thc apparent stress cr app consists of two terms: the i s defined to be a body force that constrains the velocity 

first term is a stress tensor including the external force field of the dispersions, and its explicit form is written as 



f shear = (p t (v° x (r)-v x (r))/h,0,0),r£ V. (7) 



and the second term is a stress tensor arising from the 
acceleration of thc dispersion. In experimental viscoelas- 
tic measurements of dispersions, the acceleration term can If the external force is introduced as a boundary force, 



be ignored, because the relaxation time scales related to then the development of the velocity from the bound- 
the acceleration are considerably smaller than the exper- ary to bulk (i.e., thc propagation modes) are observed 
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for short time scales, such as the kinematic time scale If we consider a particle of radius 0.4/im in water at 
T f{= PfO) /??)• In the simulation, the propagation modes room temperature, the units of space and time correspond 
influence the overall viscoelastic properties of the disper- to be 0.1pm and 0.0112 psec, respectively. In this case the 
sion. On the other hand, in experimental measurements simulated range of the frequency, / = oj/27t, is between 
the propagation modes are ignored. Using the external 7.1 and 2842 kHz. 
forces mentioned above, we can eliminate the propagation 

modes numerically. 3.1 Test of the simulation method 

To measure the storage modulus G'(u>) and loss modu- 

In order to test the validity of our method, we applied the 

lus G"(uj) 7 we monitor the xy component of the apparent 

method to a Newtonian host fluid that does not contain 

shear stress, <J^ P I P , and shear rate 7 as a function of time. 

dispersed particles. Figurc[T](a) shows the time evolutions 

In general, the xy component of the stress is written as 

of the shear rate 7 and shear stress a% pp at oj = 0.1. It 
a xy P = a o cos(wi — S) (8) can 1-^ seen ^iat both curves develop in time with the 

where CTo is the amplitude of the stress and 6 is the phase same P hase > Le ^ 5 = °- This behavior represents the typ- 
difference between &%* and 7. In our model, the shear ical features of Newtonian fluids. The loss modulus G" 



■xy 



3 Simulation results 



rate is an externally controlled parameter. By using the was then calculated using Eq. © for different frequen- 

obtaincd cto and 5, we can determine the following dy- cies - Fi § ure □ 0) displays the frequency dependence of 

namic viscoelastic moduli- ^ oss modulus G" for the host fluid. The modulus and 

„• r „_ s: frequency are non-dimensionalized by a 3 /fcsT and a 2 /Dq, 
G'(oS) = ^^uj, G"{u) = ^^u. (9) 

'° ^° respectively, where Dq is the self-diffusion coefficient of a 

Brownian particle at infinite dilution, Dq = ksT /Qirrja. 
The loss modulus increases linearly with w, and its slope 
Three-dimensional simulations were performed at several is equal to the viscosity of the host fluid. The solid line 
frequencies ranging from 0.0005 to 0.2. The amplitude of indicates the loss modulus of the host fluid, r/u. These re- 
strain, 70 = Jq^ 2u j(s)ds = 70/w, was set to be 0.2. The suits show that the viscoelastic properties of the host fluid 
spatial and temporal units are expressed in terms of the are correctly reproduced by the simulation, 
lattice spacing A and pfA 2 /rj, respectively. The system We next investigated a concentrated dispersion corn- 
lengths are L x = L y = L z = 64. The parameters of the posed of spherical particles fluctuating in the host fluid, 
simulation are a = 4, £ = 2, p p = 1, pj = 1, r\ = 1, and The volume fraction of the dispersed particles was <P = 
ksT — 7. The dispersed particles are initially randomly 0.41. Figure [2] (a) shows the time evolutions of the shear 
distributed. rate 7 and shear stress a^ P l p at cj = 0.1 and <P = 0.41. 
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The behavior of the shear stress of the dispersion differs [5] We can see the values of G' at high frequencies approach 

from that of the host fluid shown in Fig. [TJ We can see to the theoretical value. 

that there are phase differences between the shear rate t-, n , , . , j_i j_i 

r furthermore, the loss modulus is larger than the stor- 

and shear stress (i.e., S 0), and the amplitude of the , , , , , 

v ' ' " * age modulus, and the dispersion behaves like a viscous 

shear stress becomes greater than that of the shear rate. a . , TT , „ il ,. , •, , 

tluid. Here, we define the particle relaxation time to be 

From the obtained S, a and 70, we can calculate the ^ = a 2/ 6Do( $^ where Do{<p) is thc diffusion coefficient 
storage modulus G' and thc loss modulus G" at differ- of a spherical particlc in the d i S p ers ion at thermal equilib- 
ent frequencies. Figure H (b) shows the frequency depen- rimn FrQm thc cquilibrium calculations (7 = 0), we can 
dence of the storage modulus and the loss modulus for cstimatc thc diffusion coefficient at <P = 0.41, resulting in 
thc dispersion at £ = 0.41. Thc modulus and frequency ^ = r _i „ l 5 x 1Q _ 2 xhe characteristic frequency cu p is 
are scaled by a 3 /k B T and a 2 /D , respectively. For low algo indicatcd ag an arrow in Fig pjj At u „ ^ (wTp _ ^ 
frequencies, G'(u) increases linearly with and G"(w) thc onset of elastic ity is dearly observed. These behav- 
increases linearly with w. As w increases, G' grows mono- iorg accurately represent the typica l viscoelastic features 
tonically until it reaches a plateau region, while G" de- of concentrated dispcrsions m , 
vclops up to a linear region, and its slope is smaller than 

Finally, we examined the dispersion behavior of par- 

thc slope at low frequencies. The slope at high frequencies 

tides at higher concentrations, <P = 0.46 and 0.51. Thc 

represents the high frequency viscosity 77^ , which is given 

frequency dependence of G' and G" for different volume 
fractions is shown in Fig. [3J the data at <P = and 0.41 are 

?4 = lim G"(u)/u. (10) 

also plotted for comparison. For <P — 0.46, thc structures 
From Fig.[2]the high frequency viscosity was roughly esti- of the particles are still randomly distributed under shear, 
mated as 77^ ~ 5.2. Next we focus on thc high-frequency wheres for <P = 0.51 the system forms crystallize phases 
elastic shear modulus G'^, which is defined as completely. For <P = 0.46 and 0.51, each modulus shifts 

to higher values than those at <P = 0.41. At low frequen- 

Jim o G / ( W ) = G , 00 . (11) 

cies, the dependence of both moduli is remarkable: the 
We compared our results with a theoretical expression for storage modulus G' rises rapidly with increasing <P, caus- 
ey of lmrd spheres, which has been derived by Lionbcrger ing the storage modulus to become larger than the loss 
and Russel [57j. To evaluate the theoretical value, we per- modulus at low frequencies and <P = 0.51. Thc dispersion 
formed a numerical integration by using the obtained 77^ then becomes elastic. Furthermore, the loss modulus G" 
and Percus-Yevick distribution function for <P — 0.41. The at = 0.51 is a concave curve, with a minimum at low 
calculated high-frequency modulus, G'^, is shown in Fig. frequencies. 
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We found that the viscoelastic properties of the con- 
centrated dispersions depend strongly on the volume frac- 
tion. In addition, at <P = 0.51, the storage modulus in- 
creases very slowly with increasing uj, and at uj = 0.02, 
a crossover from elastic- to viscous-dominant regions is 
observed. 

3.2 Comparison of simulations and experiments 

The simulation results obtained in this work were com- 
pared with experimental results for the dynamic viscoelas- 
ticity of colloidal dispersions, which were measured by sev- 
eral group |28l29l30t3"T] . These experiments use colloidal 
dispersions of spherical particles dispersed in solvent. All 
experimental data were also plotted as a function of di- 
mensionless variables, G'a 3 /fcsT, G"a 3 /fcsT, and uja 2 /Dq. 
In FiglU the simulation results for clastic and loss modu- 
lus at <P = 0.41 were plotted together with those measured 
experimentally by Shikata and Pearson |28| . No fitting pa- 
rameters are used. The viscoelastic responses of the simu- 
lations agree well with those of the experiment, although 
the volume fraction between them is not exactly the same. 
This clearly confirms our simulation can provide the typ- 
ical viscoelastic behavior of colloidal dispersions in fluid 
states. 

Comparison of simulation results and several experi- 
ments, including concentrated dispersions with higher 
would give us a comprehensive information for the un- 
derstanding of the dynamical behavior of colloidal disper- 
sions. Figure [5] show the experimental data on the elastic 
and loss modulus of concentrated dispersions over a wide 
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range of volume fractions from 0.37 to 0.56 [28 29 30 31 , 
and our simulation results were also plotted. These ex- 
perimental data were measured over a range of too 1 / Dq 
from 10~ 5 to 10 3 , and our simulation results lie within 
the range. The experimental curve shifts at lower values of 
uia 2 /Do as the volume fraction increases and this reflects 
the characteristic time scales of system become longer, and 
the dynamics of dispersed particles slows down. 



The experimental data for high volume fraction vary 
widely depending on each experiment, and it is considered 
that the viscoelastic response becomes sensitive to the de- 
tails of direct interactions between particles. The repulsive 
potential used in this study can produce the features of 
dispersions with relatively small volume fraction (See Fig. 
[4]). The simulation at <P = 0.51 has a crystalline structure 
over the whole frequency range, where a colloidal cystal is 
formed. Thus the simulation results at $ = 0.51 are com- 
pletely different from the experiments, for which the struc- 
tures are in fluid and glasslikc states. Recently Crassous 
et al have examined the effect of crystallization on the 
viscoclasticity of the dispersions at low frequency regime 
in the vicinity of the glass transition, and they found the 
elastic modulus is greater than the loss modulus at low 
frequencies even in fluidlike states. This behavior is quali- 
tatively similar to the simulation results at <P = 0.51, and 
it represents solidlike respnses.. In the present work the 
dispersions in glassy states were not studied. 
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4 Conclusion 

We have developed a DNS method for simulating dynam- 
ics of solid particles dispersed in simple and complex fluids 
|9ll0j . This method, called SPM, has also been success- 
fully applied to simulate properties of dispersions under 
several non-equilibrium conditions, such as electrophoresis 
of charged spherical particles under external electric fields 
|19j . The method was then modified to introduce thermal 
fluctuations into the dispersions , so that one can also 
simulate situations where the thermal fluctuations and hy- 
drodynamic interactions acting among dispersed particles 
are both important. We have carried out systematic simu- 
lations for dispersions composed of spherical particles with 
and without steady shear flow, in order to analyze the dif- 
fusion process of dispersed particles in detail and to 
investigate the nonlinear viscosity of the system [55] . 

In the present paper, we report an important extension 
of the SPM for analyzing the viscoelastic properties of par- 
ticle dispersions immersed in host fluids, by introducing a 
temporally oscillatory external force into the system. To 
be consistent with the usual periodic boundary conditions, 
shear flow with a zigzag profile was employed. The valid- 
ity of the method was examined by evaluating the stor- 
age G'(u>) and loss G"(lo) moduli of a system comprising 
identical spherical particles dispersed in an incompress- 
ible Newtonian host fluid at volume fractions of # = 0, 
0.41, 0.46, and 0.51, for flow rates spanning a range of fre- 
quencies, 0.005 < u) < 0.2. We confirmed that the method 
could successfully reproduce 1) purely viscous responses 
for <P = and 2) typical viscoelastic responses for <P = 0.41 
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and 0.46, in excellent agreement with experimental data 
obtained for colloidal dispersions. 

To our knowledge, the present study is the first suc- 
cessful attempt to calculate the complex modulus of parti- 
cle dispersions using DNS-type methods. Further applica- 
tions of our DNS method to more complex systems, such 
as dispersions of non-spherical particles, particles in poly- 
mer matrices, or dispersions of aggregating particles, are 
promising. 
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Fig. 2. The complex moduli of the dispersion at $ = 0.41: (a) 
Time evolutions of the shear stress <r£y P and the shear rate 7 at 
uj — 0.1. (b) Frequency dependence of the storage modulus G' 
(O) an d loss modulus G" () for the dispersion. The modulus 
and frequency are scaled by a 3 /knT and a 2 /Do, respectively. 
G^o is the high frequency elastic modulus, which was obtained 
by the numerical integration of a theoretical expression for G'^ 
(dotted line) [27], and 77^ is the high frequency viscosity. The 
down arrow indicates the inverse of the particle relaxation time, 
u P =r- l = 6D (<P)/a 2 . 
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Fig. 3. Frequency dependence of (a) the storage modulus G' 
and (b) the loss modulus G" at different frequencies: 4> = 0.51 
(black symbols), $ = 0.46 (gray symbols), $ = 0.41 (white 
symbols), and <E> = 0.0 (*). The modulus and frequency are 
scaled by a 3 /ksT and a 2 /Do, respectively. The amplitude of 
strain is set to 0.2. 
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Fig. 4. The frequency dependent shear moduli of the simula- 
tions at <P = 0.41 (open symbol) and the experimental results 
at <S> — 0.37 for 53 nm silica particles in ethylene glycol/glycol 
(filled symbol), which were measured by Shikata and Pearson 
[28], and the high frequency elastic modulus from the simula- 
tion (dotted line). 
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Fig. 5. Data for the elastic and loss modulus of colloidal dis- 
persions in an oscillatory shear flow. The experimental data 
were collected from several papers (Shikata and Pearson [28] : 
53 nm silica sphere in ethylene glycol/glycol, Mason [33; 210 
nm silica sphere in ethylene glycol, Crassous [21] ; thermosensi- 
tive particles in water, and Watanabe [31] ; 46 nm silica sphere 
in ethylene glycol/glycol). The numbers indicate the volume 
fraction of particles. 
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